/*******************************************************************************
 * NAME
 *        mxm4 - cross product of two 4 by 4 matrices
 *
 * SYNOPSIS
 *        #include <mcce.h>
 *
 *        void mxm4(float m1[][4], float m2[][4], float m3[][4]);
 *
 * DESCRIPTION
 *        The mxm4() function calculates the cross product of matrices m1[4][4]
 *        and m2[4][4], and stores the resulted 4 by 4 matrix in m3[4][4].
 *
 * EXAMPLE
 *        #include <mcce.h>
 *
 *        int main()
 *        {  float MA[4][4] = {{ 1,  2,  3,  4},
 *                             { 5,  6,  7,  8},
 *                             { 9, 10, 11, 12},
 *                             {13, 13, 15, 16}};
 *           float MB[4][4] = {{16, 15, 14, 13},
 *                             {12, 11, 10,  9},
 *                             { 8,  7,  6,  5},
 *                             { 4,  3,  2,  1}};
 *           float MC[4][4];
 *
 *           mxm4(MA, MB, MC);
 *
 *           return 0;
 *        }
 *
 * AUTHOR
 *        Junjun Mao, 05/26/2003
 *******************************************************************************/

void mxm4(float m1[4][4], float m2[4][4], float m3[4][4])
{  m3[0][0]=m1[0][0]*m2[0][0]+m1[0][1]*m2[1][0]+m1[0][2]*m2[2][0]+m1[0][3]*m2[3][0];
   m3[0][1]=m1[0][0]*m2[0][1]+m1[0][1]*m2[1][1]+m1[0][2]*m2[2][1]+m1[0][3]*m2[3][1];
   m3[0][2]=m1[0][0]*m2[0][2]+m1[0][1]*m2[1][2]+m1[0][2]*m2[2][2]+m1[0][3]*m2[3][2];
   m3[0][3]=m1[0][0]*m2[0][3]+m1[0][1]*m2[1][3]+m1[0][2]*m2[2][3]+m1[0][3]*m2[3][3];

   m3[1][0]=m1[1][0]*m2[0][0]+m1[1][1]*m2[1][0]+m1[1][2]*m2[2][0]+m1[1][3]*m2[3][0];
   m3[1][1]=m1[1][0]*m2[0][1]+m1[1][1]*m2[1][1]+m1[1][2]*m2[2][1]+m1[1][3]*m2[3][1];
   m3[1][2]=m1[1][0]*m2[0][2]+m1[1][1]*m2[1][2]+m1[1][2]*m2[2][2]+m1[1][3]*m2[3][2];
   m3[1][3]=m1[1][0]*m2[0][3]+m1[1][1]*m2[1][3]+m1[1][2]*m2[2][3]+m1[1][3]*m2[3][3];

   m3[2][0]=m1[2][0]*m2[0][0]+m1[2][1]*m2[1][0]+m1[2][2]*m2[2][0]+m1[2][3]*m2[3][0];
   m3[2][1]=m1[2][0]*m2[0][1]+m1[2][1]*m2[1][1]+m1[2][2]*m2[2][1]+m1[2][3]*m2[3][1];
   m3[2][2]=m1[2][0]*m2[0][2]+m1[2][1]*m2[1][2]+m1[2][2]*m2[2][2]+m1[2][3]*m2[3][2];
   m3[2][3]=m1[2][0]*m2[0][3]+m1[2][1]*m2[1][3]+m1[2][2]*m2[2][3]+m1[2][3]*m2[3][3];

   m3[3][0]=m1[3][0]*m2[0][0]+m1[3][1]*m2[1][0]+m1[3][2]*m2[2][0]+m1[3][3]*m2[3][0];
   m3[3][1]=m1[3][0]*m2[0][1]+m1[3][1]*m2[1][1]+m1[3][2]*m2[2][1]+m1[3][3]*m2[3][1];
   m3[3][2]=m1[3][0]*m2[0][2]+m1[3][1]*m2[1][2]+m1[3][2]*m2[2][2]+m1[3][3]*m2[3][2];
   m3[3][3]=m1[3][0]*m2[0][3]+m1[3][1]*m2[1][3]+m1[3][2]*m2[2][3]+m1[3][3]*m2[3][3];

   return;
}
